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Summary 

The Weibull distribution has been widely adopted for the statistical description and inference of 
fatigue data. This document provides user instructions, examples, and verification for software to analyze 
gear fatigue test data. The software was developed presuming the data are adequately modeled using a 
two-parameter Weibull distribution. The calculations are based on likelihood methods, and the approach 
taken is valid for data that include type I censoring. The software was verified by reproducing results pub- 
lished by others. 


Introduction 

This document provides user instructions, examples, and verification data for software to analyze gear 
fatigue test data. The software was developed on the presumption that the data are adequately modeled 
using a two-parameter Weibull distribution. It is based on the likelihood methods described by Meeker 
and Escobar (1998). The software can be used to determine 

1. Maximum likelihood estimates of the Weibull distribution 

2. Data for contour plots of relative likelihood for two parameters 

3. Data for contour plots of joint confidence regions 

4. Data for the profile likelihood of the Weibull distribution parameters 

5. Data for the profile likelihood of any percentile of the distribution 

6. Likelihood-based confidence intervals for parameters and/or percentiles of the distribution 

The software can account for tests that are suspended without failure. The statistical terminology 
for suspended tests is “censoring.” The analytical approach for the software is valid for type I censoring, 
which is the removal of unfailed units at a prespecified time. Confidence regions and intervals are calcu- 
lated using the likelihood ratio method. Guidance for the philosophy and interpretation of statistical confi- 
dence intervals can be found in the text of Hahn and Meeker ( 1 99 1 ). 

The Weibull distribution has been widely adopted for the statistical description and inference of 
fatigue data. The Weibull cumulative distribution function F is defined as 


F(t) = 1 - exp 



where t is the time to failure, p is the shape parameter, and q is the scale parameter. For the Weibull 
distribution, the probability density function can have a variety of shapes, and the hazard function may 
be an increasing or decreasing function of time. The shapes of these curves depend only on the shape 
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table I.— data SET 1 

[Data from Krantz et al., 2000.] 


Test 

number 

Test time, 
millions cycles 

Test result 

1 

40.26 

Failure 

2 

53.94 



3 

59.88 



4 

67.68 



5 

95.76 



6 

134.22 



7 

198.48 



8 

256.20 

1 


9 

299.52 

Suspended 

10 

301.56 



11 

303.66 



12 

304.08 



13 

305.80 



14 

306.90 



15 

335.40 




TABLE II.— DATA SET 2 
[Data from Townsend and Shimski, 1994.] 


Test 

number 

Test time, 
millions cycles 

Test result 

1 

7.08 

Failure 

2 

8.15 



3 

15.52 



4 

21.34 



5 

33.37 



6 

35.70 



7 

38.31 



8 

40.26 



9 

40.74 



10 

49.47 



11 

51.70 



12 

55.29 



13 

59.17 



14 

70.18 



15 

127.07 



16 

130,95 



17 

133.38 


r 

18 

300.00 

Suspended 

19 

300.00 

Suspended 

20 

300.00 

Suspended 


parameter. Typically, surface fatigue data for gears have distributions with shape parameters less than 
three. Therefore, the probability distribution functions typically are skewed right. For a given Weibull 
distribution, one can determine the mode, mean, median, variance, and higher moments analytically by 
employing the gamma function (Cohen, 1973). 

The software developed herein was written in the FORTRAN 90 programming language. Users of the 
software are required to write short main programs that provide the data to be analyzed and also call the 
appropriate subroutines. The examples provided in this manual were compiled using the Microsoft 
PowerStation Fortran Compiler, version 4.0 (Professional Edition). 

The data analyzed as examples for this document (tables I and II) are those from gear fatigue experi- 
ments. These data have previously been analyzed using regression-based methods, and the results 
reported (Krantz et al., 2000). The data are again analyzed and reported in this document using software 
based on likelihood methods. 

Verification information for the software is provided in appendix A. The software was verified by 
reproducing results published by Meeker and Escobar (1998). Appendix B provides descriptions and 
instructions for using subroutines. Some calculations require the use of routines from a standard math 
subroutine library (IMSL, 1997; IMSL is a registered trademark of Visual Numerics, Inc.) or an equiva- 
lent. Subroutine wmaxll (appendix B) is based on the code of Keats, Lawrence, and Wang (1997). 


Estimates of Distribution Parameters 

Data Set 1 

The main program that analyzed the data of table I to provide maximum likelihood estimates of the 
Weibull distribution shape and scale parameters is presented below as source code listing 1 . Allocatable 
dimensioning is used (lines 8 and 9) to make certain that the array sizes match the value of “n,” the vari- 
able describing the number of tests (line 7). For the array “censor,” a value of “zero” is given to indicate a 
test run to failure whereas a value of “one” is given to indicate a test suspended without failure (lines 10 
to 15). A call to subroutine wmaxll provides maximum likelihood estimates of shape and scale parameters 
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(line 40). With the maximum likelihood values for the shape and scale available, subroutine wblinverse 
provides maximum likelihood estimates for any desired percentile of interest (lines 48 to 52). 

The main program produced the output presented below following source code listing 1. The maxi- 
mum likelihood estimates for shape and scale parameters are 1.037 and 376.7, respectively. 


Source Code Listing 


ine Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 


1 * main program to analyze superfinish gear fatigue test results 

2 implicit none 

3 integer n, i 

4 real, allocatable : : time(:} 

5 integer, allocatable:: censor (: ) 

6 real shape , scale , frac, life 

7 n = 15 

8 allocate (time (n) ) 

9 allocate (censor (n) ) 

10 do i=l , 8 

11 censor(i) = 0 

12 end do 

13 do i=9,n 

14 censor (i) = 1 

15 end do 

16 open (20, f ile=" superdat2 . txt" ) 

17 I HERE WE NEED TO ENTER THE ARRAY OF TEST TIMES 

18 time(l)= 40.26 

19 time (2 ) = 53.94 

20 time (3) = 59.88 

21 time (4 ) = 67.68 

22 time(5) = 95.76 ...... 

23 time (6) = 134.22 

24 time ( 7) = 198.48 

25 time (8) = 256 .20 

26 time ( 9) = 299.52 

27 time(10)= 301.56 

28 time(ll)= 303.66 

29 time (12 ) = 304.08 

30 time (13) = 305.80 

31 time (14 ) = 306.90 

32 time ( 15) = 335.40 

33 do i=l,n 

34 if (censor(i) .eq. 0 ) then 

35 write (20, 100) time(i) 

36 else if (censor(i) .eq. 1 ) then 

37 write (20, 101) time(i) 

38 end if 

39 end do 

40 call wmaxll (time , censor, n, shape , scale) 

41 write (20 , 105) 

42 write (20, 106) shape, scale 

43 100 format (f 10 . 2 , ' test failed '} 

44 101 format (f 10 . 2 , 7 test suspended ' ) 

45 105 format (' maximum likelihood estimates ') 

46 106 f ormat ( ' shape = ',el2.4, r scale = \el2.4) 

47 frac = 0.0 

48 do 13*1,9 

49 frac = frac + 0.1 

50 call wblinverse (shape , scale, frac, life) 

51 write (20, 107) frac, life 

52 end do 

53 107 format (f 6 . 2, ' proportion life = ' ,fl0.4) 

54 close (20) 

55 stop 

56 end 
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Output produced by this main program was as follows: 


40.26 

test 

failed 

53 . 94 

test 

failed 

59.88 

test 

failed 

67.68 

test 

failed 

95.76 

test 

failed 

134 .22 

test 

failed 

198.48 

test 

failed 

256.20 

test 

failed 

299.52 

test 

suspended 

301.56 

test 

suspended 

303.66 

test 

suspended 

304.08 

test 

suspended 

305 . 80 

test 

suspended 

306 . 90 

test 

suspended 

335.40 

test 

suspended 

maximum 

likelihood estimates 

shape = 

.103 

.7E+01 scale = 


10 

proportion 

life 

= 

43.0405 

20 

proportion 

life 

= 

88.7278 

30 

proportion 

life 

= 

139.4505 

40 

proportion 

life 

= 

197.1552 

50 

proportion 

life 

= 

264 . 6016 

60 

proportion 

life 

= 

346 .2899 

70 

proportion 

life 

= 

450.5648 

80 

proportion 

life 

= 

596.0463 

90 

proportion 

life 


841.8331 


Estimates of Distribution Parameters 

Data Set 2 

The main program that analyzed the data of table II to provide maximum likelihood esti mat es of the 
Weibull distribution shape and scale parameters is presented below as source code listing 2. The source 
code listing 2 mirrors source code listing 1, differing only in the lines defining the number and results of 
experiments (lines 7 to 32). The main program produced the output presented below following source 
code listing 2. The maximum likelihood estimates for shape and scale parameters are 0.8616 and 103.0, 
respectively. 


Source Code Listing 2 

Line Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 

1 I main program to analyze superfinish gear fatigue test results 

2 implicit none 

3 integer n, i 

4 real, allocatable : : time(:) 

5 integer, allocatable:: censor ( : ) 

6 real shape, scale, frac, life 

7 n = 20 

8 allocate {time (n) ) 

9 allocate (censor (n) } 

10 do i=l, 17 

11 censor (i) * 0 

12 end do 

13 do i=18 , n 

14 censor (i) = 1 

15 end do 

16 open (20, f ile="superdat 1 . txt" } 

17 I HERE WE WEED TO ENTER THE ARRAY OF TEST TIMES 

18 time (1) = 7.08 

19 time (2) = 8.15 
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20 time (3) = 15.52 

21 time (4 ) = 21.34 

22 time (5) = 33.37 

23 time (6 ) = 35.70 

24 time (7) = 38.32 

25 time (8) = 40.26 

26 time ( 9) = 40.74 

27 time (10) = 49.47 

28 time ( 11) = 51.70 

29 t ime ( 12) - 55.29 

30 time ( 13} = 59.17 

31 time { 14) = 70.81 

32 time ( 15) = 127.07 

33 time (16)= 130.95 

34 time (17) = 133.38 

35 time (18) - 300.00 

36 time (19)= 300.00 

37 time (20)= 300.00 

38 do i=l f n 

34 if (censor (i) .eq. 0 ) then 

39 write (20 , 100 ) time(i) 

40 else if (censor (i) .eq. 1 ) then 

41 write (20, 101) . time (i) 

42 end if 

43 end do 

44 call wmaxll (time , censor , n, shape , scale) 

45 write (20 , 105) 

46 write (20 , 106) shape, scale 

47 100 format (f 10 . 2 , ' test failed ') 

48 101 format (f 10 . 2 , ' test suspended ') 

49 105 format ( ' maximum likelihood estimates ') 

50 106 format (' shape = ',el2.4,' scale = \el2.4) 

51 frac =0.0 

52 do i=l , 9 

53 frac = frac + 0.1 

54 call wblinverse (shape , scale , frac, life) 

55 write (20 , 107) frac,life 

56 end do 

57 107 format (f 6 . 2 , ' proportion life = ',fl0.4) 

58 close (20) 

59 stop 

60 end 

Output produced by this main program was as follows: 

7.08 test failed 
8.15 test failed 
15.52 test failed 
21.34 test failed 
33.37 test failed 

35.70 test failed 
38.31 test failed 
40.26 test failed 
40.74 test failed 
49.47 test failed 

51.70 test failed 
55.29 test failed 
59.17 test failed 
70.81 test failed 

127.07 test failed 
130.95 test failed 
133.38 test failed 

300.00 test suspended 

300.00 test suspended 

300.00 test suspended 
maximum likelihood estimates 
shape = . 8616E+00 scale = . 1030E+03 

.10 proportion life = 7.5571 

.20 proportion life = 18.0554 

.30 proportion life = 31.1180 
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.40 proportion life 
.50 proportion life 
.60 proportion life 
.70 proportion life 
.80 proportion life 


47.2135 
67.2835 
93.0217 
127.7066 
178 . 8617 


.90 proportion life 


271.0445 


Joint Confidence Regions 

Data Set 1 

The main program that analyzed the data of table I to determine joint confidence regions for the 
simultaneous estimates of the shape parameter and 10-percent life is presented below as source code list- 
ing 3. Output files were written to provide data for likelihood ratios and confidence levels (lines 25 to 28). 
Allocatable dimensioning is used (lines 30 and 31) to make certain that the array sizes match the value of 
“n,” the variable describing the number of tests (line 29). For the array “censor,” a value of “zero” is 
given to indicate a test run to failure whereas a value of “one” is given to indicate a test suspended with- 
out failure (lines 32 to 37). The range of values for the shape and scale parameters must be provided (lines 
59 to 62). A call to subroutine wjlone provides the data required for contour plots. The subroutine writes 
results to FORTRAN units numbered 9 and 10. The output from the program is two data files, each con- 
taining 1681 lines. Each line consists of three numbers corresponding to a 10-percent life value, a shape 
parameter value, and either the likelihood ratio or the likelihood-ratio-based joint confidence value. The 
joint confidence values are displayed as a contour plot (fig. 1), providing a graphical display of plausible 
values for simultaneous estimates of the 10-percent life and the shape parameter. 



Figure 1 . — Likelihood-ratio-based joint confidence regions for data of table I. 
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Source Code Listing 3 


Line Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 


! a main program to calculate likelihood based probability contours 
! and relative likelihood contours of Weibull distributions from censored data 

i 

! this program presumes that the two parameters to 
! be plotted are the shape parameter and the 10% life 
! version 1 . 0 

! written by Tim Krantz On 2-11-2001 

implicit none ! declare all variables 

integer n I the number of tests 

real, allocatable; : time(:) ! the test times 

integer, allocatable:: censor(:) ! the censoring information 

real tenmax, tenmin I the max and min values for the 10 percent life 

! to be used in calculations 
real shapemin, shapemax ! the max and min values for the shape factor 

! to be used in calculations 


integer i I a do loop counter 

i change program as needed for a particular case 
! change only lines betweep the two lines of stars 
, *************************************************************** 

! the outputs for this program will be written to these two files 
i the file formats are as follows {each has 3 columns) 

\ column 1 is the shape parameter 
I column 2 is the 10 percent life 

I column 3 is the likelihood ratio or the confidence number as a percentage 
l file for unit #10 is the confidence number data 
l file for unit #9 is the likelihood ratio data 
OPEN (10, FILE = w super -prob. txt") 

OPEN (9, FILE = "super-ratio.txt") 
n = 15 

allocate (time (n) ) 
allocate (censor (n) ) 
do i=l , 8 
censor (i) = 0 
end do 


do i = 9 , n 
censor (i) = 1 


end do 

I HERE WE NEED TO ENTER THE ARRAY OF TEST TIMES 
time (1) = 40.2600 
time (2)= 53.9400 
t ime ( 3 ) = 59.8800 
t ime (4 ) = 67.6800 
t ime ( 5 ) = 95.7600 
time (6) = 134.220 
time (7) = 198.480 
time (8) = 256.200 
time (9)= 299.520 
time ( 10) = 301.560 
time (11)= 303.660 
time (12 ) = 304.080 
time (13)= 305.800 
time(14)= 306.900 
time ( 15) = 335.400 


l here we set the bounds for 10 percent life and shape factors 
I tenmax is the largest 10 percent life value 
I tenmin is the smallest 10 percent life value 
I shapemax is the largest shape value 
I shapemin is the smallest shape value 


tenmax = 140. 
tenmin - 2.0 
shapemax = 2 . 
shapemin = 0.4 

I none of the lines below need to be changed to run a particular case 

j ******************************************************************** 

! the call to this subroutine checks for errors in the censoring array 


call ccheck (n, censor) 

I now we ask for the solution 

call wj lone (time, censor, n, shapemax, shapemin, tenmax, tenmin) 
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Data Set 2 


69 close(lO) 

70 close (9) 

71 stop 

72 end 


The main program that analyzed the data of table II to determine joint confidence regions for the 
simultaneous estimates of the shape parameter and 10-percent life is presented below as source code list- 
ing 4. The source code listing 4 mirrors source code listing 3. The output from the program is two data 
files, each containing 1681 lines. Each line consists of three numbers corresponding to a 10-percent life 
value, a shape parameter value, and either the likelihood ratio or the likelihood-ratio-based joint confi- 
dence value. The joint confidence values are displayed as a contour plot (fig. 2), providing plausible 
values for simultaneous estimates of the 1 0-percent life and the shape parameter. 



Figure 2. — Likelihood-ratio-based joint confidence regions for data of table II. 


Source Code Listing 4 

Line Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 

1 ! a main program to calculate likelihood based probability contours 

2 ! and relative likelihood contours of Weibull distributions from censored data 

3 i 

4 I this program presumes that the two parameters to 

5 i be plotted are the shape parameter and the 10% life 
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6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 


! version 1.0 

! written by Tim Krantz On 2-11-2001 

implicit none l declare all variables 

integer n ! the number of tests 

real, allocatable : : time ( : ) ! the test times 

integer, allocatable:: censor (: ) 1 the censoring information 

real tenmax, tenmin I the max and min values for the 10 percent life 

! to be used in calculations 
real shapemin, shapemax ! the max and min values for the shape factor 
! to be used in calculations 


integer i ! a do loop counter 

' change program as needed for a particular case 
! change only lines between the two lines of stars 

i *************************************************************** 

i the outputs for this program will be written to these two files 
I the file formats are as follows {each has 3 columns) 

1 column 1 is the shape parameter 
I column 2 is the 10 percent life 

l column 3 is the likelihood ratio or the confidence number as a percentage 
I file for unit #10 is the confidence number data 
I file for unit #9 is the likelihood ratio data 
OPEN (10, FILE * "base-prob.txt") 

OPEN (9, FILE = "base-ratio.txt") 
n = 20 


allocate (time (n) ) 
allocate (censor (n) ) 
do i=l , 17 
censor (i) = 0 
end do 


do i=18,n 
censor (i) = 1 
end do 

I HERE WE NEED TO ENTER THE ARRAY OF TEST TIMES 
time(l)= 7.08 
time (2 ) * 8.15 
time (3 ) = 15.52 
time(4)= 21.34 
time (5) = 33.37 
time (6) = 35.70 
time (7) = 38.32 
time ( 8 ) = 40.26 
time{9)= 40.74 
time (10) = 49.47 
time (11) = 51.70 
time (12 ) = 55.29 
time (13 ) = 59.17 
time (14 ) = 70.81 
time (15) = 127.07 
time (16)= 130.95 
time (17) = 133.38 
time (18)= 300.00 
time (19) = 300.00 
time (20) = 300.00 

here we set the bounds for 10 percent life and shape factors 
tenmax is the largest 10 percent life value 
tenmin is the smallest 10 percent life value 
shapemax is the largest shape value 
shapemin is the smallest shape value 


tenmax = 30. 
tenmin = 1.5 
shapemax = 1.4 
shapemin = 0.4 

I none of the lines below need to be changed to run a particular case 

1 ******************************************************************** 

I the call to this subroutine checks for errors in the censoring array 
call ccheck (n, censor) 

! now we ask for the solution 

call wj lone (time , censor, n, shapemax, shapemin, tenmax, tenmin) 
close (10) 
close (9) 
stop 


77 end 
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Profile Likelihood and Profile-Likelihood-Based Confidence Intervals 


Data Set 1 

The main program that determined the profile likelihood and profile-likelihood-based confidence 
intervals for the 10-percent life for the data of table I is presented below as source code listing 5. The 
main program determines the profile likelihood over a range of values for a specified percentile of inter- 
est, in this case the 10th percentile. Any percentile can be investigated by changing one line in the code 
for variable “pfrac” (line 49). Allocatable dimensioning is used (lines 18 and 19) to make certain that the 
array sizes match the value of “n,” the variable describing the number of tests (line 17). For the array 
“censor,” a value of “zero” is given to indicate a test run to failure whereas a value of “one” is given to 
indicate a test suspended without failure (lines 24 to 29). Array “time” (lines 31 to 45) provides the test 
times. The subroutine plvpercent is called (line 56) to determine the relative likelihood ratio for a speci- 
fied 10-percent life. By calling the subroutine repeatedly within a do-loop (lines 50 to 59), data for a 
profile-likelihood plot are calculated. The output from the main program was plotted using two vertical 
scales to show both the profile-likelihood ratio values and the confidence level (fig. 3). Table III lists 
some profile-likelihood ratios and the corresponding confidence level. A 90-percent confidence interval 
for the 10-percent life is shown graphically on figure 3. The lower and upper bounds can be determined 
more precisely by running the program again with an appropriate range of values for the 10-percent life 
(figs. 4 and 5). 


TABLE III. — CONFIDENCE LEVELS 
AND CORRESPONDING PROFILE- 
LIKELIHOOD RATIOS 


Confidence level 

Ratio 

0.50 

0.7965 

.60 

.7018 

.70 

.5844 

.80 

.4399 

.85 

.3548 

.90 

.2585 

.95 

.1465 

.97 

.0811 

.99 

.0362 



Figure 3. — Profile-likelihood ratios and confidence levels for 1 0-percent life for data of 
table I. 


NASA7TM — 2002-21 1 1 09 


10 





7 8 9 10 11 


1 0-percent life, millions of cycles 

Figure 4. — Profile-likelihood-based confidence interval for 10-percent life for data of 
table I (lower bound). 



10-percent life, millions of cycles 


Figure 5.— Profile-likelihood-based confidence interval for 10-percent life for data of 
table I (upper bound). 


The subroutine plvpercent called by the main program determines the shape and scale factors corre- 
sponding to the largest likelihood value by an iterative search method. This method requires that maxi- 
mum and minimum values for the shape parameter be provided to begin the search. The minimum and 
maximum values used in this example were hardcoded in the subroutine as 0.05 and 8.0 respectively. 
Users who wish to analyze data sets that require inspection of shape values outside the just-stated range 
need to modify the subroutine (appendix B). 


Source Code Listing 


Line Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 


1 ! a main program to calculate likelihood based profile 

2 I likelihood curves of Weibull distributions from censored data 
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3 » 

4 l this program presumes that a percentile 

5 l parameter is the one of interest 

6 ! written by Tim Krantz On 2-13-2001 

7 implicit none I declare all variables 

8 I the following variables represent inputs 

9 l they are hard-coded in this main program 

10 integer n ! the number of tests 

11 real, allocatable r : time ( : ) I the test times 

12 integer, allocatable:: censor {: ) ! the censoring information 

13 integer i < do loop counter 

14 real blv, spblv, scblv 

15 real pfrac 

16 real tmin,tmax,dt 

17 n = 15 

18 allocate (time (n) ) 

19 allocate (censor (n) ) 

20 I HERE WE NEED TO SUPPLY CENSORING INFORMATION 

21 1 A VALUE OF "0" MEANS THAT A FAILURE OCCURRED 

22 ! A VALUE OF "1" MEANS THAT TEST WAS SUSPENDED (CENSORED) WITH NO FAILURE 

23 ! ENTER THE VALUES AS INTEGERS 

24 do i*l,8 

25 censor (i) = 0 

26 end do 

27 do i=9,n 

28 censor(i) = l 

29 end do 


30 

I HERE WE 

NEED TO ENTER THE ARRAY OF TEST TIMES 

31 

time ( 1) = 

40 .2600 

32 

time (2 ) = 

53 . 9400 

33 

time (3) = 

59 . 8800 

34 

time (4 ) = 

67.6800 

35 

time (5) = 

95 . 7600 

36 

time (6 ) = 

134.220 

37 

time (7) = 

198 .480 

38 

time ( 8 ) = 

256 .200 

39 

time ( 9) = 

299 . 520 

40 

time (10) = 

301 . 560 

41 

time ( 11 ) = 

303 . 660 

42 

time ( 12 ) = 

304 . 080 

43 

time (13) = 

305.800 

44 

time (14) = 

306 . 900 

45 

time ( 15) * 

335.400 


46 l the call to this subroutine checks for errors in the censoring array 

47 call ccheck2 (n, censor) 

48 open (10 , FILE= "super_L10_profile.txt") 

49 pfrac = 0.1 

50 tmin = 4. 

51 tmax = 135. 

52 dt = (tmax-tmin) / 100. 

53 tmin = tmin - dt 

54 do i=l , 101 

55 tmin = tmin + dt 

56 call plvpercent (time , censor , n, tmin, pfrac , blv, spblv, scblv) 

57 write (10 , 103) tmin, spblv, scblv, blv 

58 103 format (4el6 . 7) 

59 end do 

60 stop 

61 end 


Data Set 2 

The main program that determined the profile likelihood and profile-likelihood-based confidence 
intervals for the 10-percent life for the data of table II is presented below as source code listing 6. The 
main program mirrors that of source code 5. The output from the main program was plotted using 
two vertical scales to show both the profile-likelihood ratio values and the confidence level (fig. 6). 

A 90-percent confidence interval for the 1 0-percent life is shown graphically on figure 6. 
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Figure 6.— Profile-likelihood ratios and confidence levels for 10-percent life for data of 
table II. 


The subroutine plvpercent called by the main program determines the shape and scale factors corre- 
sponding to the largest likelihood value by an iterative search method. This method requires that maxi- 
mum and minimum values for the shape parameter be provided to begin the search. The minimum and 
maximum values used in this example were hardcoded in the subroutine as 0.05 and 8.0 respectively. 
Users who wish to analyze data sets that require inspection of shape values outside the just- stated range 
need to modify the subroutine (appendix B). 

Source Code Listing 6 

Line Source Line Microsoft Fortran PowerStation Compiler. Version 4.0 

1 I a main program to calculate likelihood based profile 

2 l likelihood curves of Weibull distributions 

3 I from censored data 

4 : 

5 l this program presumes that a percentile parameter is the one of interest 

6 l version 1 . 0 

7 l written by Tim Krantz On 2-13-2001 

8 implicit none ! declare all variables 

9 ! the following variables represent inputs 

10 I they are hard-coded in this main program 

11 integer n ! the number of tests 

12 real, allocatable : : time ( : ) I the test times 

13 integer, allocatable : : censor (:) ! the censoring information 

14 integer i ! do loop counter 

15 real blv, spblv, scblv 

16 real pfrac 

17 real tmin,tmax,dt 

18 n = 20 

19 allocate (time (n) ) 

20 allocate (censor (n) ) 

21 ! HERE WE NEED TO SUPPLY CENSORING INFORMATION 

22 i A VALUE OF "0" MEANS THAT A FAILURE OCCURRED 

23 ! A VALUE OF "1" MEANS THAT TEST WAS SUSPENDED (CENSORED) WITH NO FAILURE 

24 do i=l,17 

25 censor (i) = 0 

26 end do 

27 do i=18 , n 

28 censor (i) = 1 

29 end do 

30 I HERE WE NEED TO ENTER THE ARRAY OF TEST TIMES 
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31 

t ime { 1 ) = 

7.08 

32 

time (2) = 

8.15 

33 

time (3) = 

15.52 

34 

time (4 ) = 

21.34 

35 

time (5) = 

33,37 

36 

time (6) = 

35.70 

37 

time (7) = 

38.32 

38 

time (8) = 

40.26 

39 

time (9) = 

40 . 74 

40 

time (10) = 

49.47 

41 

time (11)— 

51.70 

42 

time ( 12 ) = 

55,2 9 

43 

time (13) = 

59.17 

44 

time (14) * 

70.81 

45 

time (15) = 

127.07 

46 

time (16) - 

130 . 95 

47 

time (17) = 

133.38 

48 

time (18) = 

300 . 00 

49 

time (19) = 

300.00 

50 

time (20) = 

300 . 00 


51 \ the call to this subroutine checks for errors in the censoring array 

52 call ccheck2 (n, censor) 

53 open (10, FILE^ "basel_L10_profile.txt") 

54 pfrac - 0.10 

55 write (6,*) ' enter min value ' 

56 read (5,*) tmin 

57 write ( 6 ,*) ' enter max value ' 

58 read (5,*) tmax 

59 dt = (tmax-tmin) / 100. 

60 tmin = tmin - dt 

61 do i=l,101 

62 tmin = tmin + dt 

63 call plvpercent (time , censor , n ( tmin, pfrac, blv # spblv, scblv) 

64 write {10 , 103 ) tmin, spblv, scblv, blv 

65 103 format (4el6. 7} 

66 end do 

67 stop 

68 end 


Weibull Plot With Confidence Intervals 

Profile-likelihood-based confidence interv als can be determined for any percentile of interest using 
main programs similar to source code listings 5 and 6. The variable “pfrac” determines the percentile to 
be investigated. For example, using “pfrac = 0.30” one can calculate data for the 30-percent life. Collec- 
tions of such confidence intervals were calculated, and the results are provided in figures 7 and 8. These 
two figures are plotted using Weibull coordinates so that the Weibull cumulative distribution frequency 
will plot as a straight line. The test data points are plotted at the positions of exact median ranks (Jaquelin, 
1993). 
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Figure 7.— Weibull plot for data of table I with profile-likelihood-based confidence 
intervals. Intervals shown are pointwlse 90-percent confidence intervals of percentile 
estimates. 



Cycles to failure, millions 

Figure 8. — Weibull plot for data of table II with profile-likelihood-based confidence 
intervals. Intervals shown are pointwise 90-percent confidence intervals of percentile 
estimates. 
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Summary of Results 


Software has been developed to analyze gear fatigue test data. The source code listings discussed 
provided examples illustrating how to define the test data to be analyzed and how to call the appropriate 
subroutines to analyze the data. The software can be used to determine 

1 . Maximum likelihood estimates of the Weibull distribution 

2. Data for contour plots of relative likelihood for two parameters 

3. Data for contour plots of joint confidence regions 

4. Data for the profile likelihood of the Weibull distribution parameters 

5. Data for the profile likelihood of any percentile of the distribution 

6. Profile-likelihood-based confidence intervals for parameters and/or percentiles of the distribution 

Contour plots of relative likelihood (not illustrated in this document) can be obtained from the output 
files as provided by source code listings 3 and 4. Profile-likelihood calculations for Weibull distribution 
parameters were not illustrated directly. A point for the profile-likelihood plot of the shape parameter can 
be determined using subroutine wl ratio (appendix B). A profile likelihood for the scale parameter can be 
obtained using a source code in the manner of listings 5 and 6 and analyzing for the 63.2 percentile 
(“pfrac” = 0.632). 
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Appendix A 

Verification of Software 


The software was verified by reproducing results published by Meeker and Escobar (1998). The main 
programs to produce the figures in this appendix (figs. 9 to 14) were similar to the source code listings 
found in the main body of this document. 



Figure 9. — Contours of equal likelihood ratios. Weibull shape parameter p = 1/a 
and Weibull scale parameter r\ = exp (pi). In the manner of figure 8.1 from Meeker 
and Escobar (1 998). 



Figure 10. — Contours of equal joint confidence levels. Weibull shape parameter 
P = 1/a and Weibull scale parameter n] = exp (pi). In the manner of figure 8.4 
from Meeker and Escobar (1 998). 
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Profile likelihood 



9.8 10.0 10.2 10.4 10.6 10.8 


Figure 1 1 .-Profile-likelihood method to determine confidence interval for scale 
parameter. Weibull scale parameter r| = exp (|x). In the manner of figure 8.5 from 
Meeker and Escobar (1998). 



a> 
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Figure 12. — Profile-likelihood method to determine confidence interval for shape 
parameter. WeibuIT shape parameter P =T/a. In the manner of figure 876 from 
Meeker and Escobar (1998). 
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Confidence level 



6000 8000 10 000 12 000 14 000 16 000 18 000 20 000 

1 0-percent life, fo.l 

Figure 13. — Contours of equal joint confidence levels for Weibull shape parameter 
and 10-percent life, fo.l - Weibull shape parameter p = 1/a. In the manner of 
figure 8.7 from Meeker and Escobar (1 998). 



6000 8000 10 000 12 000 14 000 16 000 18 000 20 000 


10-percent life, fo.l 

Figure 14.— Profile-likelihood method to determine confidence interval for 1 0-percent 
life, #o.i- In the manner of figure 8.8 from Meeker and Escobar (1998). 
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Appendix B 

Subroutine Descriptions and User Instructions 

Subroutine CCHECK 

A subroutine to do “error” checking on censoring array, all values should be either 1 or 0 integers. 

Usage 

subroutine ccheck (n, censor) 

Arguments 

n (integer) size of array “censor” (input) 

censor (integer array) array of dimension [n], censoring information (input) 

Comments 

If all elements in array “censor” equal either 0 or 1, no action is taken. 

If any element does not meet the criteria, a line of text is written to FORTRAN units numbered 6, 9, 
and 10, and then execution of the main program is halted. 

Other subroutines called 

None 
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Subroutine CCHECK2 


A subroutine to do error checking on censoring array. All values should be either 1 or 0 integers. 

Usage 

subroutine ccheckl (n, censor) 

Arguments 

n (integer) size of array “censor’ (input) 

censor (integer array) array of dimension [n], censoring information (input) 

Comments 

If all elements in array “censor” equal either 0 or 1 , no action is taken. 

If any element does not meet the criteria, a line of text is written to FORTRAN units numbered 6 and 
10, and then execution of the main program is halted. 

Other subroutines called 

None 
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Subroutine FIXSHAPEMAXLL 


A subroutine to calculate the maximum likelihood estimates of the scale parameter for a two- 
parameter Weibull distribution, for a presumed shape and censored data. 

Usage 

subroutine fixshapemaxll (times, censor, n, shape, scale) 


Arguments 


times 

n 

censor 

shape 

scale 


(real array) 
(integer) 
(integer array) 
(real) 

(real) 


array of dimension [n], test times (input) 
size of array “censor” (input) 

array of dimension [n], censoring information (input) 
presumed Weibull shape parameter (input) 
calculated Weibull scale parameter (output) 


Comments 

This routine is based on the code published by Keats and Lawrence (1997). Although arguments 
are passed as single precision, calculations within the subroutine are done in double precision. 

For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

Other subroutines called 


None 
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Subroutine PLVPERCENT 


A subroutine to calculate a point on a profile-likelihood curve for a Weibull distribution at a particu- 
lar cumulative distribution function (CDF) percentile of interest and a presumed value for the time to 
failure at that percentile. This routine can be called repeatedly using different values for the assumed 
time to failure (atime) to generate data for a profile-likelihood curve. 

Usage 


subroutine phpercent (time, censor, n, atime, pfrac, ratio, spblv, scblv) 

Arguments 


time 

(real array) 

array of dimension [n], test times (input) 

censor 

(integer array) 

array of dimension [n], censoring information (input) 

n 

(integer) 

size of array “censor” (input) 

atime 

(real) 

assumed time to failure (input) 

pfrac 

(real) 

assumed percentile of CDF corresponding to 
value of atime (input) 

ratio 

(real) 

profile ratio for best likelihood value [blv] (output) 

spblv 

(real) 

shape factor corresponding to blv (output) 

scblv 

Comments 

(real) 

scale factor corresponding to blv (output) 


For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

Data are rescaled w ithin the subroutine to avoid extreme values. 

The Weibull shape parameter resulting in the profile-likelihood value is found within the subroutine 
by an iterative method. The shape value is found to within a tolerance of 0.01. A range for the shape must 
be provided to start the iterative process, and the shape parameter resulting in the profile-likelihood value 
must be contained within that range. The subroutine uses a range from 0.05 to 8.0. These values could be 
adjusted to fit the needs of a particular data set, but one might encounter numerical problems for extreme 
scale parameter values since it is used as an exponent in calculations. The routine writes a line of text to 
unit number 6 if the routine fails to converge within 1000 iterations. 

Other subroutines called 


wmaxll 

wlratio 
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Subroutine WBLINVERSE 


A subroutine to calculate an inverse Weibull cumulative distribution function (CDF) value. 


Usage 
subroutine u 

? bl inverse (shape, scale, frac, life) 

Arguments 

shape 

(real) 

Weibull shape parameter (input) 

scale 

(real) 

Weibull scale parameter (input) 

frac 

(real) 

Weibull CDF percentile of interest (input) 

life 

(real) 

time to failure for provided percentile of CDF 

Comments 

None 

Other subroutines called 


None 
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Subroutine WFLSHAPE 


A subroutine to calculate Weibull distribution likelihood-based profile data (likelihood ratio). 

Usage 

subroutine wflshape (z, censor, n, shapemax, shapemin) 

Arguments 

z 

censor 
n 

shapemax 
shapemin 

Comments 

For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test 
suspended without failure (type I censoring). 

No outputs are returned to the calling program. However, output is written to FORTRAN unit 
number 9. If a file for unit number 9 was not yet opened, the user will be prompted for a filename. 

The likelihood profile ratio is calculated for 41 values of the shape factor, the values equally spaced 
between the values passed as shapemax and shapemin. The output written to unit number 9 consists of 
two columns of numbers of fixed length, where the numbers are the shape factor and the calculated likeli- 
hood ratio. 

Data are rescaled within the routine to avoid extreme values. 

Other subroutines called 

wmaxll 

fixshapemaxll 
w l ratio 


(real array) 
(integer array) 
(integer) 

(real) 

(real) 


array of dimension [n], test times (input) 

array of dimension [n], censoring information (input) 

size of array “censor” (input) 

maximum shape value of range to be calculated (input) 
minimum shape value of range to be calculated (input) 
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Subroutine WJLONE 


A subroutine to calculate Weibull distribution likelihood-based joint confidence ratios and corre- 
sponding confidence levels to be plotted as contour plots. This routine presumes that the Weibull distribu- 
tion is parameterized by the shape parameter and 10-percent life. 

Usage 

subroutine wjlone (z, censor, n, shapemax, shapemin, tenmax, tenmin) 


Arguments 

array of dimension [n], test times (input) 
array of dimension [n], censoring information (input) 
size of array “censor” (input) 

maximum shape value of range to be calculated (input) 
minimum shape value of range to be calculated (input) 
maximum 10-percent life to be calculated (input) 
minimum 10-percent life to be calculated (input) 

Comments 


z 

censor 

n 

shapemax 

shapemin 

tenmax 

tenmin 


(real array) 
(integer array) 
(integer) 

(real) 

(real) 

(real) 

(real) 


For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

No outputs are returned to the calling program. However, output is written to FORTRAN units 
numbered 9 and 10. If files were not yet opened, the user will be prompted for filenames. 

The likelihood profile ratios and confidence levels are calculated for a full array of 41x41 values 
(1681 lines) of the shape parameter and 10-percent life, the values equally spaced between the ranges 
defined by the values passed as shapemax to shapemin and tenmax to tenmin, respectively. The output 
written to unit number 9 consists of three columns of numbers of fixed length, the numbers being the 
shape parameter, the 10-percent life, and the calculated likelihood ratio. The output written to unit number 
10 consists of three columns of numbers of fixed length, where the numbers are the shape parameter, the 
10-percent life, and the calculated joint confidence level. 

Data are rescaled within the routine to avoid extreme values. 


NOTE: Confidence levels less than 10 percent are written to the output file as equal to 10 percent. 
Confidence levels greater than 99.99 percent are written to the output file as equal to 99.99 percent. 

Other subroutines called 


wmaxll 

wlratio 

chidf (IMSL routine — 


see IMSL (1997); IMSL is a registered trademark of Visual Numerics, Inc.) 


NASA/TM— 2002-21 1 109 


27 



Subroutine WJLTWO 


A subroutine to calculate Weibull distribution likelihood-based joint confidence ratios and corresponding 
confidence levels to be plotted as contour plots. This routine presumes that the Weibull distribution is 
parameterized by the shape parameter and scale parameter (or 63.2-percent life). 

Usage 

subroutine wjltwo (z, censor, n, shapemax, shapemin, scalemax, scalemin) 

Arguments 

z ’ 

censor 
n 

shapemax 
shapemin 
scalemax 
scalemin 

Comments 

For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

No outputs are returned to the calling program. However, output is written to FORTRAN units 
numbered 9 and 10. If Fdes were not yet opened, the user will be prompted for filenames. 

The likelihood profile ratios and confidence levels are calculated for a full array of 41x41 values 
(1681 lines) of the shape parameter and scale parameter, the values equally spaced between the ranges 
defined by the values passed as shapemax to shapemin and scalemax to scalemin, respectively. The out- 
put written to unit number 9 consists of three columns of numbers of fixed length, the numbers being the 
shape parameter, the scale parameter, and the calculated likelihood ratio. The output written to unit num- 
ber 10 consists of three columns of numbers of fixed length, where the numbers are the shape parameter, 
the scale parameter, and the calculated joint confidence level. 

Data are rescaled within the routine to avoid extreme values. 

NOTE: Confidence levels less than 10 percent are written to the output file as equal to 10 percent. 
Confidence levels greater than 99.99 percent are written to the output file as equal to 99.99 percent. 

Other subroutines called 

wmaxll 

wlratio 

chidf (IMSL routine — see IMSL (1997); IMSL is a registered trademark of Visual Numerics, Inc.) 


(real array) 
(integer array) 
(integer) 

(real) 

(real) 

(real) 

(real) 


array of dimension [n], test times (input) 

array of dimension [n], censoring information (input) 

size of array “censor” (input) 

maximum shape value of range to be calculated (input) 
minimum shape value of range to be calculated (input) 
maximum scale value of range to be calculated (input) 
minimum scale value of range to be calculated (input) 
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Subroutine WLOGLV 


A subroutine to calculate the log-likelihood value of a single life test. Can be called repeatedly, with 
the outputs summed, to calculate the joint log-likelihood value of a data set. 

Usage 

subroutine wloglv (sc, sh, t, fail, llv) 


Arguments 


SC 

(real) 

Weibull scale parameter (input) 

sh 

(real) 

Weibull shape parameter (input) 

t 

(real) 

test time (input) 

fail 

(integer ) 

flag indicating test to failure of test censored (input) 

llv 

(real) 

log-likelihood value (output) 


Comments 

For variable fail, a value of 0 indicates a test to failure whereas a value of 1 indicates a test suspended 
without failure (type I censoring). 

Other subroutines called 

None 
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Subroutine WLRATIO 


A subroutine to calculate the numerator of a likelihood ratio of a data set using presumed Weibull 
shape and scale parameters. 

Usage 

subroutine \v! ratio (time, censor, n, shape, scale, lv) 

Arguments 

time 
censor 
n 

shape 
scale 
lv 

Comments 

For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

Other subroutines called 

wloglv 


(real array) 
(integer array) 
(integer) 

(real) 

(real) 

(real) 


array of dimension [n], test times (input) 

array of dimension [n], censoring information (input) 

size of array “censor” (input) 

presumed Weibull shape parameter (input) 

presumed Weibull scale parameter (input) 

likelihood value, or value of numerator (output) 
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subroutine WMAXLL 


A subroutine to calculate the maximum likelihood estimates of the scale parameter for a two- 
parameter Weibull distribution, for a presumed shape and censored data. 

Usage 

subroutine wmaxll (times, censor, n, shape, scale) 

Arguments 

times 
n 

censor 
shape 
scale 

Comments 

This routine is based on the code published by Keats, Lawrence, and Wang (1997). Although argu- 
ments are passed as single precision, calculations within the subroutine are done in double precision. 

For array “censor,” a value of 0 indicates a test to failure whereas a value of 1 indicates a test sus- 
pended without failure (type I censoring). 

The routine finds the maximum likelihood value for the shape parameter iteratively, where the shape 
parameter is the root of an equation. The routine is considered as converged if the step size of the 
Newton-Raphson method is less than 0.0003. If the routine fails to converge after 100 iterations, a line is 
written to FORTRAN unit number 15. 

Other subroutines called 

None 


(real array) array of dimension [n], test times (input) 

(integer) size of array “censor” (input) 

(integer array) array of dimension [n], censoring information (input) 

(real) calculated Weibull shape parameter (output) 

(real) calculated Weibull scale parameter (output) 
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